      subroutine vumat(
C Read only (unmodifiable)variables -
     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
     2  stepTime, totalTime, dt, cmname, coordMp, charLength,
     3  props, density, strainInc, relSpinInc,
     4  tempOld, stretchOld, defgradOld, fieldOld,
     5  stressOld, stateOld, enerInternOld, enerInelasOld,
     6  tempNew, stretchNew, defgradNew, fieldNew,
C Write only (modifiable) variables -
     7  stressNew, stateNew, enerInternNew, enerInelasNew )
C
      include 'vaba_param.inc'
C
      dimension props(nprops), density(nblock), coordMp(nblock,*),
     1  charLength(nblock), strainInc(nblock,ndir+nshr),
     2  relSpinInc(nblock,nshr), tempOld(nblock),
     3  stretchOld(nblock,ndir+nshr),
     4  defgradOld(nblock,ndir+nshr+nshr),
     5  fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
     6  stateOld(nblock,nstatev), enerInternOld(nblock),
     7  enerInelasOld(nblock), tempNew(nblock),
     8  stretchNew(nblock,ndir+nshr),
     8  defgradNew(nblock,ndir+nshr+nshr),
     9  fieldNew(nblock,nfieldv),
     1  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     2  enerInternNew(nblock), enerInelasNew(nblock)
C
      character*80 cmname
C
      parameter ( zero = 0.d0, one = 1.d0, two = 2.d0,three=3.d0,four=4.d0,five=5.d0,six=6.d0,d_thresh = 0.99d0,
     *     third = 1.d0 / 3.d0, half = 0.5d0)

      parameter (nblkLocal = 136)

      real*8 STRESSM,PSTRESS(6)
      real*8 SEQ,SEQ1

      real*8 SYIELD,SYIELD_eps,SYIELD_Temp
      real*8 SYIELD_T

      real*8 RHS

      real*8 EIrate           !应变率
      real*8 EI               !等效塑性应变
      real*8 DEI,DEI1              !等效塑性应变增量
      real*8 DEImin,DEImax  

      real*8 dlamda

      real*8 E,V,SYIELD0,H,ppow,C_jc,G,Lamda,mpow,Tm,T0,alphaT
      real*8 de_expan
      common/EI/EI(19809)
      common/DEI/DEI(19809)

      E = PROPS(1)
      V = PROPS(2)

      SYIELD0 = PROPS(3)

      H = PROPS(4)

      ppow = PROPS(5)

      C_jc = PROPS(6)

      mpow = PROPS(7)

      Tm = PROPS(8)

      T0 = PROPS(9)

      alphaT = PROPS(10)




      do k = 1,nblock
        n_ele = INT(stateNew(k,10))  
        SYIELD_T = stateOld(k,22)
C         EI = stateOld(k,1)
C         DEI = stateOld(k,2)
        EIrate = DEI(n_ele)/dt

        G = E/two/(one+V)
        Lamda = V*E/(one+V)/(one-two*V)

        de_expan = alphaT*(tempNew(k) - tempOld(k))

        trace = strainInc(k,1)+strainInc(k,2)+strainInc(k,3) - three*de_expan
        stressNew(k,1) = stressOld(k,1) + Lamda*trace + two*G*(strainInc(k,1) - de_expan)
        stressNew(k,2) = stressOld(k,2) + Lamda*trace + two*G*(strainInc(k,2) - de_expan)
        stressNew(k,3) = stressOld(k,3) + Lamda*trace + two*G*(strainInc(k,3) - de_expan)
        stressNew(k,4) = stressOld(k,4)               + two*G*strainInc(k,4)
        stressNew(k,5) = stressOld(k,5)               + two*G*strainInc(k,5)
        stressNew(k,6) = stressOld(k,6)               + two*G*strainInc(k,6)

C     静水压力
        STRESSM = (stressNew(k,1)+stressNew(k,2)+stressNew(k,3))/three
C     偏应力
        PSTRESS(1) = stressNew(k,1) - STRESSM
        PSTRESS(2) = stressNew(k,2) - STRESSM
        PSTRESS(3) = stressNew(k,3) - STRESSM
        PSTRESS(4) = stressNew(k,4)
        PSTRESS(5) = stressNew(k,5)
        PSTRESS(6) = stressNew(k,6)
C     试应力等效应力
        SEQ = SQRT(((stressNew(k,1)-stressNew(k,2))**two+(stressNew(k,1)-stressNew(k,3))**two
     1      +(stressNew(k,2)-stressNew(k,3))**two+six*(stressNew(k,4)**two+stressNew(k,5)**two
     2      +stressNew(k,6)**two))/two)

C     屈服应力
        if ( EIrate .lt. one ) then
              SYIELD_eps = one
        else 
              SYIELD_eps = one + C_jc*log(EIrate)
        end if

        if ( tempNew(k) .lt. T0 ) then
            SYIELD_Temp = one
        else if (tempNew(k) .gt. Tm) then 
            SYIELD_Temp = zero
        else 
            SYIELD_Temp = one - ((tempNew(k) - T0)/(Tm - T0))**mpow
        end if

        SYIELD = (SYIELD0 + H*EI(n_ele)**ppow)*SYIELD_eps*SYIELD_Temp
        if ( SYIELD .lt. SYIELD_T )SYIELD=SYIELD_T 

        if ( SEQ .gt. SYIELD ) then
              RHS = SEQ - SYIELD 

              nn = zero

              DEImin = zero
              DEImax = SEQ/three/G
              DEI1 = zero
              do while(abs(RHS) .gt. 1e-8)

                nn=nn+1

                DEI1 = half*(DEImin+DEImax)

                EIrate = DEI1/dt 

                if ( EIrate .lt. one ) then
                      SYIELD_eps = one
                else 
                      SYIELD_eps = one + C_jc*log(EIrate)
                end if

                SYIELD = (SYIELD0 + H*(EI(n_ele) + DEI1)**ppow)*SYIELD_eps*SYIELD_Temp
                if ( SYIELD .lt. SYIELD_T )SYIELD=SYIELD_T 

                RHS = SEQ - three*G*DEI1 - SYIELD

                if ( RHS .gt. zero ) DEImin = DEI1
                if ( RHS .lt. zero ) DEImax = DEI1
                if ( RHS .eq. zero ) goto 100

                if ( nn .gt. 100 ) then
                    write(6,*)"too much iterations"
                    !call XPLB_EXIT
                    goto 100
                end if

              end do 
100           continue  

            EI(n_ele) = EI(n_ele) + DEI1

            DEI(n_ele) = DEI1
            stateNew(k,1) = EI(n_ele) 
            stateNew(k,2) = DEI(n_ele) 


            if ( SEQ .gt. zero ) then
                dlamda = three/two*DEI(n_ele)/SEQ
            else 
                dlamda = zero
            end if

            stressNew(k,1) = stressNew(k,1) - (Lamda*(PSTRESS(1)+PSTRESS(2)+PSTRESS(3))+two*G*PSTRESS(1))*dlamda
            stressNew(k,2) = stressNew(k,2) - (Lamda*(PSTRESS(1)+PSTRESS(2)+PSTRESS(3))+two*G*PSTRESS(2))*dlamda
            stressNew(k,3) = stressNew(k,3) - (Lamda*(PSTRESS(1)+PSTRESS(2)+PSTRESS(3))+two*G*PSTRESS(3))*dlamda
            stressNew(k,4) = stressNew(k,4) - two*G*PSTRESS(4)*dlamda
            stressNew(k,5) = stressNew(k,5) - two*G*PSTRESS(5)*dlamda
            stressNew(k,6) = stressNew(k,6) - two*G*PSTRESS(6)*dlamda
            


        end if
        stateNew(k,22)=SYIELD

C Update the specific internal energy -
        stressPower = half * (
     1    ( stressOld(k,1)+stressNew(k,1) )*strainInc(k,1)
     1    +     ( stressOld(k,2)+stressNew(k,2) )*strainInc(k,2)
     1    +     ( stressOld(k,3)+stressNew(k,3) )*strainInc(k,3)
     1    + two*( stressOld(k,4)+stressNew(k,4) )*strainInc(k,4) 
     1    + two*( stressOld(k,5)+stressNew(k,5) )*strainInc(k,5)
     1    + two*( stressOld(k,6)+stressNew(k,6) )*strainInc(k,6))
C
        enerInternNew(k) = enerInternOld(k)
     1    + stressPower / density(k)
C
C Update the dissipated inelastic specific energy -
      SEQ1 = SQRT(((stressNew(k,1)-stressNew(k,2))**two+(stressNew(k,1)-stressNew(k,3))**two
     1      +(stressNew(k,2)-stressNew(k,3))**two+six*(stressNew(k,4)**two+stressNew(k,5)**two
     2      +stressNew(k,6)**two))/two)

      plasticWorkInc = SEQ1*DEI(n_ele)

        enerInelasNew(k) = enerInelasOld(k)
     1    + plasticWorkInc / density(k)

        stateNew(k,15)=enerInelasNew(k)



      end do  

      return
      end

      subroutine vusdfld(
c Read only variables -
     1   nblock, nstatev, nfieldv, nprops, ndir, nshr, 
     2   jElem, kIntPt, kLayer, kSecPt, 
     3   stepTime, totalTime, dt, cmname, 
     4   coordMp, direct, T, charLength, props, 
     5   stateOld, 
c Write only variables -
     6   stateNew, field )
c
      include 'vaba_param.inc'
c
      dimension jElem(nblock), coordMp(nblock,*), 
     1          direct(nblock,3,3), T(nblock,3,3), 
     2          charLength(nblock), props(nprops), 
     3          stateOld(nblock,nstatev), 
     4          stateNew(nblock,nstatev),
     5          field(nblock,nfieldv)
      character*80 cmname
c
c     Local arrays from vgetvrm are dimensioned to 
c     maximum block size (maxblk)
c
      parameter( nrData=6 )
      character*3 cData(maxblk*nrData)
      dimension stress_prin(maxblk*nrData), jData(maxblk*nrData)
c


       do kblock = 1, nblock
            stateNew(kblock,10)=jElem(kblock)

       end do


       
c
      return
      end